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ABSTRACT 

Most analytic work to date on protostellar disks has focused on those in isolation from their environ- 
ments. However, observations are now beginning to probe the earliest, most embedded phases of star 
formation, during which disks are rapidly accreting from their parent cores and cannot be modeled in 
isolation. We present a simple, one-zone model of protostellar accretion disks with high mass infall 
rates. Our model combines a self-consistent calculation of disk temperatures with an approximate 
treatment of angular momentum transport via two mechanisms. We use this model to survey the 
properties of protostellar disks across a wide range of stellar masses and evolutionary times, and make 
predictions for disks' masses, sizes, spiral structure, and fragmentation that will be directly testable 
by future large-scale surveys of deeply embedded disks. We define a dimensionless accretion-rotation 
parameter which, in conjunction with the disk's temperature, controls the disk evolution. We track 
the dominant mode of angular momentum transport, and demonstrate that for stars with final masses 
greater than roughly one solar mass, gravitational instabilities are the most important mechanism as 
most of the mass accumulates. We predict that binary formation through disk fission, fragmentation 
of the disk into small objects, and spiral arm strength all increase in importance to higher stellar 
masses. 

Subject headings: accretion disks — binaries: general - stars: formation — -ISM xlouds 



1. INTRODUCTION 

A young star system's visible T Tauri or Herbig stage 
is preceded by a deeply enshrouded phase of rapid ac- 
cretion in which its character multiplicity, disk prop- 
erties, and tendency to form planets - is first forged. 
Although this embedded phase is likely the one during 
which most accretion onto the star occurs, the properties 
of disks during this period have received relatively little 
attention. This phase is difficult to model analytically 
because embedded disks are subject to large perturba- 
tions in the form of rapid accretion of mass and angular 
momentum, making local models and stability analyses 
problematic. Due to the high obscuration characteris- 
tic of this phase, disks are accessible primarily via ra- 
dio and submillimeter observations, and such techniques 
provide limited sensitivity and angular resolution com- 
pared to what can be achieved for T Tauri and Herbig 
AE star disks using shorter wavelengths. Our knowledge 
of massive protostellar disks is particularly limited by 
this problem, since they do not have a significant unob- 
scured phase, probably due to the destructive effects of 
ionizing radiation. These difficulties are compounded by 
the fact that massive stars form more rarely, and there- 
fore tend to lie farther away. Detections of rotation and 
infall in a few systems hint at the presence of disks during 
the embed ded phase, but are only preliminary (see re cent 
reviews by|Cesaroni et al.|2006| [20071 |Beuther|2006| and 
Beuther et al.||2007p . 

While our knowledge of the embedded phase today is 
limited, it will soon come into sharp focus as new instru- 
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ments such as the Expanded Very Large Array (EVLA) 
and Atacama Large Millimeter Array (ALMA) become 
operational. In order to predict what these telescopes 
will discover about the formation of stars across a very 
broad mass range, ~ 1 — 12OM0, we present evolutionary 
models of star-disk systems reacting to infall at very high 
rates. We concentrate our efforts on the physical pro- 
cesses that control disk evolution, such as the torque from 
a turbulent infall, the reprocessing of starlight by the in- 
fall envelope, and the character of the self-gravitational 
instabilities. The disk itself we model with a highly 
simplified, single-zone treatment. Although multidimen- 
sional simulations provide a much more detailed view of 
disk f ormation and evolution during the embedded phase 
(e.g. |Goodwin et al.||2004| |Krumholz et al.||2007b|), their 



high computational cost and the limited range of physical 
processes they include mean that these simulations can 
explore only small regions of parameter space. They can- 
not easily make predictions across a broad range of stel- 
lar mass scales and evolutionary times. In contrast, our 
semi-analytic approach permits us to incorporate more 
physical effects and explore the consequences of environ- 
mental parameters more rapidly, and in a more system- 
atic way. This serves two complementary ends: the the- 
oretical goal of understanding angular momentum trans- 
port and fragmentation in the embedded phase, and the 
observational goal of making concrete predictions about 
the properties of young, massive disks. 

Although we include a range of physical processes in 
our models; the most important in driving the evolution 
of embedded disks in our calculations is self-gravity. Self- 
gravity plays a central role in mediating angular momen- 
tum transport and triggering fragmentation into a binary 
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or multiple system. Its importance in star formation has 



2.1. 



1984), and our study is 



ations which incor porate 



long been recognized ( Larson 
preceded by evolutionary calcu 

accretion and self-gravity into one-d imensio nal ( Lin fc 
|Pringle|JT987l |l990l |Nakamoto fc Nakagawa||1994[ |l995l 
Hueso & Guillot 2005[) and two-dimensional ( jVorobyov 
fc ljasu||2005n2lX)6f ^imulations. Although lower m res- 
olution^ our approach is distinguished from these works 
in several ways: 



i. In contrast to all one-dimensional calculations to 
date, we account for the dependence of gravita- 
tional torques on the disk-to-total mass ratio in 
addition to Toomre's instability parameter. 

ii. We consider the possibility that disks will fragment 
if they become sufficiently unstable. 

Hi. We consider fluctuations of the vector angular mo- 
mentum in the infall due to realistic turbulence in 
the collapsing cloud core. 

iv. We employ a realistic model for the irradiation of 
the disk midplane, in which starlight is reprocessed 
at the inner wall of an outflow cavity while inflow 
is occurring. 

v. We survey the conditions of intermediate-mass and 
massive star formation, rather than focusing exclu- 
sively on conditions in nearby low-mass star form- 
ing regions such as Taurus. 

The first of these is important for all protostellar disks, 
since the disk mass is never negligible when the Toomre 
parameter is small. Features (ii) and (Hi) are of primary 
importance in the formation of massive stars, which ac- 
crete from strongly turbulent regi ons ( |Myers fc Fuller 
|1992| [McLaughlin fc Pudritz||1 997l and wEIch are likely 
to und 



lergo disk fragmentation ( |Kratter fc Matzner|2006 



hereafter KM06). Irradiation is most importan t for low" 
mass stars, w hose disks it strongly stabilizes (Matzner 
fc Levin|2005[), but it remains significant in massive star 



formation as well. 

In S|2]we outline our model for the infall rate of matter 
and angular momentum. We develop a model for disk 
accretion and fragmentation in Sj3] In f]4] we define the 
key environmental variables that control protostellar disk 
evolution and sketch a qualitative evolutionary sequence 
based on their fiducial values. In fj5]we present the results 
for our fiducial case, and explore the effect of varying our 
parameters. In ^6] we discuss the observable predictions 
that our model makes, and finally in § [7] we summarize 
our main results. 

2. INFALL ONTO DISKS 

Since we are interested in the behavior of a disk that is 
subject to strong perturbations from its environment, we 
begin building our model by constructing a prescription 
for the infall of matter and angular momentum onto a 
disk. This accretion comes from a background "core" 



(Shu 



1977] |McLaughlin fc Pudritz|p97] |McKee fc Tan" 



2003), whose properties and interaction with a disk we 



discuss in this section. 



Star Formation By Core Collapse 

We model the accretion of mass and angular momen- 
tum using the two-component core model of McKee & 
Tan| ( |2003[ ), which is a generaliz ation of the TNT (ther - 
mal plus non-thermal) model of Myers & Fuller (19921. 
In this model, which we summarize here for convenience, 
the density distribution within a core follows a two- 
component power law distribution 
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where R c is the outer radius of the core, c s c is the ther- 
mal sound speed within it (assumed to be consta nt), and 
p s is the dens ity at the core's surface. We follow |McKee 



& Tan (20031 in adopting k p — 1.5 as the fiducial value 
ot the turbulence-supported density index. Physically, 
the first term describes an envelope supported primar- 
ily by turbulent motions, while the second describes a 
thermally- supported region at its center. A model of this 
sort is fully specified in terms of the three parameters: 
the core mass 
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surface density 
and temperature 
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(3) 
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where m = 3.9 x 10 -24 g is the mean particle mass in a gas 
of molecular hydrogen and helium mixed in the standard 
cosmic abundance. Observed regions of star formation 
contain cores with masses ~ 1 — 100 M , surface densities 
S c w 0.03 — 1 g cm" 2 , and temperatures of 10 to 50 K 



(Johnstone et al. 



2001] |Plume et al.(l997 l 



The core is taken to be in approximate hydrostatic 
balance initially, and this condition specifies the required 
level of turbulent support. The non-thermal velocity dis- 
persion in the shell at radius r is 



a(rf = 
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GM(r) 
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where M(r) is the mass at radii of r or less and 4>b — 2.8 
approximately accounts for the magnetic contribution to 
the total pressure. Except when M(r) <§C 1 M or E c < 
0.1 g cm -2 , the first term is much larger than the second, 
so that the velocity dispersion is primarily non-thermal. 

Core collapse commences at time zero, and a mass shell 
initially at radius r falls onto the disk in a time compa- 
rable to the free fall time evaluated at its initial density, 



ts(r) = [37r/3 2G'p(r)1 1 / 2 . In practice, we use the McKee 
& Tan ( 2003| ) accretion rate approximation to determine 
M on to the star-disk system as a function of the to- 
tal core mass and the current amount of mass that has 
accreted: 
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where M*r is the final disk plus stellar mass, M*d is the 
current disk plus stellar mass, tg. s is the free-fall time 
evaluated at the core surface (i.e. at p — p(R c )), 



= 3(2 - k p ) 
q 2(3 - k p ) ' 



(7) 
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E C:0 = S c /(gcm 2 ), and </>*, </>* n th, and 0* t h are con- 
stants of order unity that depend on the polytropic index 
and magnetic field strength. The efficiency factor 



M r 



(9) 



represents the fraction of the core mass that lands on the 
star-disk system rather than bein g blown out by proto- 
stellar outflows. We again follow |McKee fc Tan| ( 2003 1 
in adoptin g e = 0.5, a value typic al of low-mass star 
formation ( TMatzner fc McKee||2000 1. 



2.2. Angular Momentum of the Infalling Material 

Equation (I6| gives the mass infall rate M ln {f) from 
the core as afunction of time. The second component 
of our core model is to specify the corresponding rate 
of angular momentum infall 3 m (t). We compute this in 
several steps. First, we approximate the vector specific 
angular momentum j(r) averaged over a shell of mate- 
rial at radius r as described below. Then we compute 
M(t) = f^M[ n (t')/edt', the total mass from the core 
that has cither fallen onto the star-disk system or been 
ejected at time t. From the core density profile (equa- 
tion [lj) we also compute M(r) = f£ inr' 2 p(r') dr' , the 
mass of the initial core enclosed within radius r. Assum- 
ing the core accretes inside-out, we set M(r) = M(t) 
and solve for r(t), which gives the initial radius r of the 
shell of mass that arrives at the star-disk system at time 
t. Assuming that the specific angular momentum of the 
gas does not change before it reaches the disk, the an- 
gular momentum accretion rate is then simply given by 
jjn(i) = Min(i)j(r(i))._ 

The remaining step is to specify how we estimate j(r). 
Star forming cores are often modeled as solid body rota- 
tors characterized by the ratio j3 of rotational to gravi- 
tational energy, but we adopt a more realistic model in 
which t urbulent fluctuations affec t the infalling g as. Fol - 
lowing [Burkert fc B odcnhe imer| ( |2000| ), |Fisher| Q2004| ), 
Matzner fc Levin| ( 2005[ ) and KM06, we assume that 



the o bserved angular momenta of cores (Goodman et al. 
1993 1 can be modeled using an idealized turbulent ve - 
locity field. Using the method of Dubinski et al. |1995 1, 
we generate a numerical realization of an isotropic, ho- 
mogeneous, Gaussian random velocity field v(r). We re- 
quire that the power spectrum P{k) of this turbulent 
field reproduce the scalings required by turbulent sup- 
port against gravity: <r(r) 2 cx GM(r)/r cx r 2_,£p at large 
radii, so that a(r) cx r 1 / 4 for k p = 3/2. Parceval's theo- 
rem or dimensional analysis then require P(k) cx fc~ 3 / 2 . 

We note that numerical simulations of supersonic tur- 
bulence consistentl y show the steeper spectral index —2 
( Porter et al.|1992| which is understood as the spectrum 
of an individual shock and as the exact limit of Burgers 



turbulence. |Matzner| ( |2007[ ) and |Nakamura fc Li| ( |2007 1 
have shown that a shallower index is expected when tur- 
bulence is driven by protostellar outflows, however, and 
our chosen power spectrum is consistent with the line 
width-size relation for massi ve cores (e.g. 



ers 1995; Plume et al. 



1997J). 



|Caselli fc My- 
Our homogeneous velocity 
field is surely an idealization, but not a grave one. 

After scaling our numerical domain to match the core 
radius R c , we normalize v such that the one-dimensional 
velocity dispersion of a spherical shell with this radius 
equals er(r) defined in equation ([5]). In practice we fit 
R c within a 256 3 section of a 1024 3 grid of velocities, 
because periodicity causes artifacts on scales larger than 
about 1/4 of the box size. From this field we calcu- 
late the specific angular momentum j = r x v at every 
point and the mean specific angular momentum j(r) in a 
shell at radius r. Note that KM06 calculate the expected 
magnitude and dispersion of j(r) for velocity fields of pre- 
cisely this type; our results agree with their predictions 
to about 50%, which is within the scatter they predict. 

3. DYNAMICS OF THE DISK 
3.1. Approach to disk evolution 

Given the rate at which mass and angular momentum 
accrete, we must calculate the reaction of the disk. At 
any given time, our star-plus-disk system is characterized 
by the disk mass Md, the central star mass M„, and the 
total angular momentum content of the disk 3d- Given 
these quantities, and the rates of mass and angular mo- 
mentum infall Mj n and Ji n , we wish to compute the time 
rate of changes Md, M*, and 3d- 

Using the separation between the thermal, orbital, and 
accretion timescales, we assume our disks are in a ther- 
mal steady state and draining at a rate determined by 
their current global parameters. We shall later refer to 
this as the assumption of thermal and mechanical equi- 
librium. 

In § |3.2[ we estimate the disk accretion rate onto the 
central star due to various angular momentum trans- 
port mechanisms. In § |3.3| we discuss thermal equilib- 
rium in the disk, which together with the aforementioned 
condition of mechanical equilibrium allows us to self- 
consistently compute the accretion rate from the disk 
to the star M*. In § 

Finally, in 



we describe the corres pond ing 
angular momentum evolution 3d- 
§ |3.6| we discuss our prescriptions for disk fragmentation 
and binary formation. 

It is helpful before proceeding to define two dimcn- 
sionless parameters that characterize the strength of the 
disk's self-gravity. These are the disk-to-total mass ratio 

M d 



M d + M» 



(10) 



(11) 



and Toomre s (1964) instability parameter 

C S K 

where k is the radial epicyclic frequency, c s speed of den- 
sity waves, and is the disk's mass surface density. In 
practice we evaluate Q using k — * ft — {GM tot / R^) 1 ' 2 , 
the total orbital frequency, since the difference between 
them is only marginally significant even when the disk 
mass is quite large. Here Rd — i^/GM tot . We also 
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3M so that 
and il at the 



approximate c s using the isothermal sound speed. To 
characterize gravitational instability, we use the mini- 
mum value of Q - the value at Rd, the outer bound- 
ary of our active disk. In this evaluation we a ssume a 
profile oc r _1 (a choice we justify in 
^d(Rd) — Md/i^Rd)' an d we evaluate c s 
edge of the disk. 

We base our models on the fundamental assumption 
that the self-gravitational behavior of an accretion disk 
depends primarily on the structural parameters \x and 
Q - so th at it s evolution is controlled by heating and 
cooling (§ 3.3 I, which alter Q, and accretion onto and 
through the disk, which alters /x. This approach permits 
us to treat the disk's mechanical and thermal properties 
separately, before combining them into a model for its 
evolution. This division also guides our use of published 
work, since it implies that simulations with adiabatic 
equations of state and those with an imposed cooling 
rate may be combined into a mechanical model for disk 
evolution, which we may then use to model irradiated 
protostellar disks. Finally, it prompts us to treat the on- 
set of disk fragmentation and disk fission as boundaries 
in the space of fi and Q, rather than in terms of a crit- 
ical cooling rate (which is the natural and conventional 
description for si mulat ions that include cooling but not 
irradiation). In § 3.5 we argue that these descriptions 



are effectively equivalent 

Models based on this assumption are guaranteed to be 
somewhat approximate, because a disk's mechanical evo- 
lution must, at some level, reflect additional parameters: 
its dimensionality, its equation of state, the specifics of 
its heating and cooling processes, and the magnitude of 
external perturbations (like tides), to name a few. How- 
ever we expect our results to be valid, both because we 
believe that [i and Q are indeed the most significant pa- 
rameters for gravitational instability, and because our 
model is calibrated to realistic numerical models. Addi- 
tional simulations will be required to test this approach. 

3.2. Angular Momentum Transport and Disk Accretion 

A key element of our model is a prescription for angu- 
lar momentum transport and the rate M* at which mat- 
ter accretes onto the central star - or more specifically, 
the dimensionless rate M*/(Md£l). In practice we first 
construct a model for an effective [Shakura & Sunyaev 
a parameter, which we define through the steady-state 
relation 

M. - §f (12, 



so that 



M, 
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where the factor of of 3/8 comes from the assumption 
that the disk surface density falls as r _1 . We do not mean 
to imply by this that trans port induced by gravitationa l 
instability is purely local ( |Balbus fc Papaloizou 19991, 
although this does appear to be the case for sufficientl y 
thin and light disks ( |Gammie|20"o"i~l |Lodato fc Rice|2005 l. 

We divide a into two contributions: omri, due to the 
magnetorotational instability (MRI), and ckgIi due to 
gravitational instability. In keeping with the strategy 
described in § |3.1[ we consider aci to be a pure function 



of fi and Q. We combine it and the MRI contribution 
linearly: 

a = a G i + a MRI . (14) 

We create our model for agi(Q, n) using resu l ts from the 
simulations of Laughlin fc Rozyczka| fl!996||, |Rice et al. 



jgOOiyLodato fc Rice| ( |2004] , |Lodato fc Rice| ( |2005p and 

We a 



GaTrmiieT POOl). 
as discussed below. 



adopt a constant value for cxmri 



3.2.1. Overview of Simulations 

The three sets of simulations span a large fraction of 
our param eter space in Q and \x. The global disk sim- 
ulations of Laughlin & Rozyczka ( 1996 1 explore Q > 1 
and non-negligible values of /j, using a two-dimensional 
hydrodynamic, self-gravity code; they supress local frag- 
mentation by impo sing an adiabat ic equation of state. 
The simulations of Gammie (2001) represent the Hmit 
fi — > 0, for values of Q which approach unity from above, 
and are most directly applicable to quasar disks. |Gam-| 
|mie| imposes cooling with a fixed cooling time, r c , which 
is proportional to the orbital time. He finds a regime 
of steady gravity-induced turbulence, for disks that cool 
over many orbits. If r c is too short (< 30" 1 ), however, 
the disk fragments as Q drops below unity. The disk vis- 
cosity is highest at the boundary of fragmentation. An- 
gular momentum transport in this regime is quite local, 
with an effective value of a that is inversely proportional 
to the cooling rate. Our t hird source is the global SPH 
simul ations of Rice et a l. (2003), Lodato & Rice (200' 
and Lodato fc Rice|p005| , in which a cooling time cxTT 
is imposed locally; these cover the entire parameter space 
in fi. In these simulations Q is initially 2, but it descends 
towards unity. Here again, the disk fragments if Qt c is 
too small, althoug h the c ritical value of this parameter 
is different than IGammiel found. 

3.2.2. Accretion Model 

To derive a relatively simple analytic fit to the simula- 
tion data, we must extract a characteristic olq\, Q, and /i 
from the simulations listed above. Because the numeri- 
cal approaches are varied, we are unable to use the same 
method for each. The values are derived as follows for 
each type of simulation: 



1. We estimate aci from Laughlin & Rozyczka ( 1996 1 
using their equations 24-26; Q and fj, are given. 



2. Values of a from Rice et al. ( 2003 ), Lodato & 
Rice| ( |2004[ ) , I Lodato fc Rice| ( |2005| )~are taken di- 
rectly from plots, when available. Because a varies 
with radius, we take an approximate value from the 
outer region of their disk before the density begins 
to fall off steeply. When plots are not available, 
we use the critical value of r c to calculate aci at 
the fragmen tatio n boundary, which we take to be 
Q = 1 (see 



3.51. Again, fj, is given. 



3. 



Gammie 
for a dis 



: (2001 ) provides one value of «gi at Q = 1 
k with n — > 0, which we adopt. 

These values of aci axe shown in Figure [l] according to 
the estimated values of /i and Q that accompany each 
of them. We treat them as a data set to be fit within 
our analytical model for ogi, which is displayed in the 
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underlying contours in that figure. Imposing the realistic 
condition that aci is continuous and equals zero for Q > 
2 (when the gravi t ationa l instability should shut off as 
suggested by Griv (2006 1 ), we find that two components 
are required: 
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cohort = max 
and 



.).:i4[^-l)(i- / M ll \o 



aiong = max 
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(15) 
(16) 
(17) 



fact we apply equation (15 1 only to the region 



In 

Q > 1. Because we expect the gravitational torque to 
saturate when fragmentation occurs, we hold a constant, 
for a given fi, when Q < 1; this amounts to replac- 
ing Q — > max(Q, 1) in the above equations. This has 
no practical consequences for our calc ulati ons, however, 
since our treatment of fragmentation (§ 
disks from sampling values of Q much 



3.5 1 prevents our 
jelow unity. 

Our nomenclature in equation ( 15 1 reflects our inter- 
pretation. The "short" component a s hort dominates for 
Q < 1.3, hence for relatively thin disks. We think of it as 
arising from modes with rela tively high spatial wavenum - 
bers and short wavelengths ( Lodato fc Rice[2004 2005). 
Note that i ts fu nctional form resembles the model of ILinl 
& Pringle ( 1990 [ l (their eq. 16) modified by a mild /i de- 
pendence, which is comparable to the scale height depen - 



dence derived in equation (2.5) of Lin & Pringle (19871 
The "long" component ai ong is important m thicker 
disks whose instabilities are likely to be dominated by 
loosely wound, m = 2 spiral patterns. We require it be- 
cause we include t he adiabatic simulations of Laughlin 



& Rozyczka ( 1996 ) , which sample higher values of Q bc- 
cause they cannot cool. (Indeed, Q rises du ring these 



simulations.) Our fundamental assumption (§ |3.1[ ) leads 
us to incorporate these results into a single model for 
a(/i, Q), despite the difference in thermal physics. Future 
simulations can test this assumption by imposing heating 
(via irradiation, say) as well as cooling: our model im- 
plies that the derived aci will be comparable to [Laughlin 
|fc Rozyczkaf s, when Q and \i tak e simila r val ues. While 
simula tions such as Boley et al.| ( |2006[ ) and Cai et al. 



([2007 
mode 



are making dramatic process towards accurately 
ling heating, cooling, and irradiation, a wider pa- 



rameter space is neces s ary fo r com parison. We note that 



Sellwood & Carlberg 



Griv (20061 also find 



- ((1984|) and 

non-axisymmetric instabilities for massive disks with Q 
in the range 1.3 — 2. 

As shown in Figure [l] our model for aoi agrees reason- 
ably well with data from the simulations, though we fail 
to fit a couple points at ver y high jj, and low Q. No t e that 
aci for these points from |Laughlin fc Rozyczka ( 1996 1 



are uncertain themselves 

It is important to bear in mind that our accretion 
model is only a rough representation of the numerical 
results on which it is based, and that it can be improved 
as more simulations become available. For instance, we 
place no stock in the weak divergence of ai ong as /i — > 0: 
this feature is a product of our fit to numerical results at 
larger [i, and it would be an unwise extrapolation to use 




Fig. 1. — Contours of the visco sity parameter log(oQi) due to 
gravitational instabilities (eq. |15) ; white squares are contour la- 
bels. Results from numerical simulations are marked with circles, 
diamonds, and tria ngles. Circles show simulatio ns with adiabatic 
equations of state (iLaughlin & Rozyczka 1996), diamonds show 
simulations with an impose d cooling rate that reach steady state 
jLodato Sc R ice 2004 2005), and triangles show the maximum ojqi 
achieved m simulatio ns with imposed cooling that probe the frag- 
mentation boun dary QGa mmic | |2001| |Rice et a,l.]|2003| |Lodato fc| 
|Rice|2004| |2005| . Note t hat t he point at ii = corresponds to the 
purely local simulation of|Gammic ( 200lJ . The Q = 1 boundary is 
marked with a dashed line. 

our model for disks with very low /i and moderate Q. It 
does not affect our results, as our disks do not sample 
this regime. 

Finally, we assume the disk is sufficiently ionized to 
support magnetic turbulence, and we represent the MRI 
with the constant value ajyiRi = 10~ 2 . The typical v alue 
of q;mri is rather uncertain; see |Pessah et al. J2007|) for 
a synthesis of recent work, and [Hu^soTrTfuTiIot| ( |2005| ) 
for a recent consideration of observational constraints in 
low-mass protostellar disks. Gravitational torques ex- 
ceed those from the MRI for much of the accretion phase. 
We discuss the influence of «mri on our results in § |5.5| 

Figure [2] illustrates our model for the dimensionless 
accretion rate M*/ \M,0) as a function of Q and \i. We 
draw attention to several key features of the plot. First, 
note that at low Q there is a tongue-like feature that 
increases in intensity with increasing disk mass. This is 
due to the strong dependence of cohort on both Q and [i 
At higher values of Q and lower values of /x the contours 
steepen due to the weak divergence of aq on g as /i — > 0, 
which is probably not physical. The curvature towards 
higher Q and fi shows the dominance of the MRI for 
Q > 2. The fact that the dimensionless accretion rate 
takes numerical values up to 10 -2 4 , with a typical values 
~ 10 -3 ' 5 , implies that massive disks drain on timescales 
ranging from ~ 40 to a few thousand orbits, with five 
hundred orbits being typical. 

3.3. Disk Thermal Equilibrium 

We have now specified the rate at which the disk ac- 
cretes onto a central star as a function of Q and \x. How- 
ever, this does not fully specify the accretion rate, be- 
cause while /i may be directly computed from our "prim- 
itive" variables M4, M*, and J 4, the Toomre stability pa- 
rameter Q cannot be, because it depends on the sound 
speed c s and thus the temperature within the disk. We 
can determine this by requiring that the disk be in ther- 
mal equilibrium. 
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Fig. 2. — Contours of the dimensionless accretion rate from the 
disk onto the star (M* /M^Q) from all transport components of our 
model. The lowest contour level is 10 -4 ' 8 , and subsequent contours 
increase by 0.3 dex. The effect of each transport mechanism is ap- 
parent in the curvature of the contours. At Q < 1.3 the horizontal 
"tongue" outlines the region in which short wavelength instability 
dominates accretion. The more vertical slope of the contours at 
lower fi and Q > 1.3 shows the dominance of the long wavelength 
instability. The MRI causes a mild kink in the contours across the 
Q = 2 boundary and is more dominant at hi ghe r disk masses due 
to our assumption of a constant a: equation [13] illustrates that a 
constant a will cause higher accretion rates at higher values of fi. 

To compute the disk's thermal state, we follow the ap- 
proach of KM06, in which disks are heated by a combina- 
tion of stellar irradiation and viscous dissipation due to 
accretion. In equilibrium, the disk midplane temperature 
satisfies the approximate relation 



aT 4 = 



(18) 



where F v is the rate of viscious dissipation per unit area 
in the disk, F; rr is the flux of starlight (whether direct or 
reprocessed) onto the disk surface, and tr,p = kr^pH/2 
are the Rosseland and Planck optical depths to the mid- 
plane. The viscous dissipation per unit area is 



urn 



(19) 



and w e compute the opacities using the Semenov et al.| 
(20031 model for Kn y p(T): in particular, we use their 



homogeneous-aggregate dust model with normal sili- 
cates, calculated at the appropriate density. 

Low-mass stars' luminosities are accretion-dominated 
in the main accretion phase, but those above about 15 
Mq reach the zero-age main sequence (ZAMS) while 
still accreting. To include both accretion luminosity and 
other sources in our calc ulation of Jj rr , we use the proto - 
stellar evolut ion code of IKrumholz fc Thompson| ( |2007[ | , 
based on the |McKee fc Tan| ( |2003[ ) protostellar evolution 
model, to compute the luminosity of the central star 
as a function of its accretion history. The model includes 
contributions to the protostellar luminosity from accre- 
tion on the stellar surface, Kelvin-Hclmholtz contraction, 
and, once the temperature rises high enough, deuterium 
and then hydrogen burning. 

During the infall, dust within the infall envelope re- 
processes starlight and casts it down on the disk. By 
performing ray tracing within an inflow envelop e with 
a central outflow cavity, Matzner & Levin ( 20051 deter- 
mine the fraction of light received by the disk assuming 



the infall envelope is optically thick to the stellar radia- 
tion, and optically thin to its own IR re-radiation: they 
find 

p _ 0- 1 L * m] 

The weak dependence on the accretion efficiency e arises 
from a picture in which the outflow clears a fraction 
1 — e of the core, so that infall streamlines originate 
from regions separated fr om the axis by a n gles 9 such 
that cos 9 > e. Recently, Rodriguez et al. (20051 have 
observed an outflow near an U-type protostar with an 
opening angle of approximately 25°; this is in reasonable 
agreement with the model chosen here, since infall occurs 
at wider angles. 

Once the core has accreted entirely and the envelope 
can no longer re-process starlight, we make an (unreal- 
istically) abrupt switch to a model in which M m = 0. 
The star continues to acquire mass from the disk, which 
represents a non-negligible reservoir. From this point on 
we cal culate Fj rr in the manner of Chiang & Goldreich 
(19971. We first identify the fraction of L* intercepted 



by the surface which is optically thick to stellar pho- 
tons, assuming for this purpose that H oc R 9 ^ 7 and that 
the dust density is a Gaussian, of scale height H, in the 
height above the midplane. We also calculate the equi- 
librium temperature of dust in this reprocessing layer. 
We then calculate Fi rr as that fraction of the reprocessed 
radiation which is reabsorbed by the disk, allowing for 
the possibility that the disk will be optically thin at the 
relevant wavelengths. We find the reprocessing height is 
slightly larger than a scale height (1.5 H being typical); 
higher values are typical of more massive disks, which 
are more opaque. 

Though negligible during the accretion phase, we also 
include a background radiation field due to the cloud 
(modeled as an optically thin dust layer) and the cosmic 
microwave background. This prevents disks from becom- 
ing unrealistically cold at large radii and late times. Our 
cloud irradiation serves as a stand-in for one neglected 
heat source in clusters: irradiation from surrounding 
stars. This effect is likely important for (a) very dense 
regions, and (b) late times when disk radii stretch out 
to 10 4 AU. Due to the wide variance in the strength of 
this effect, we do not address heating by neighbors here. 
There is also minor heating due to the accretion shock 
that feeds the disk; however, KM06 have argued that this 
is negligible in general. 

While our background heating is only important at 
late times, we do not report results for t >2 Myr as this 
may exceed the lifetime of gas d isks, even the low-mass 
ones ( Jayawardhana et al. 20061. The uncertainties in 
our procedure therefore have little effect on the results 
we obtain. 

We have now fully specified the conditions of thermal 
and mechanical equilibrium for this disk, and we can use 
them to compute the accretion rate. Equations ( 13 1 and 
(18 1 constitute two equations for the unknowns Q and 
M*. For any given M^, A/,, and Jd, we may solve them 
to determine M*. This in turn also specifies the rate of 
change of the disk mass 

M d = M- m -M,. (21) 
Note that Md can also be modified by disk fragmentation 
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and binary formation, as described in § |3.5|and § 3.6 



3.4. An Outer Disk and the Braking Torque 

When describing standard steady-state disks, one im- 
plicitly assumes that when angular momentum is trans- 
ported radially, it travels out to large radii in an insignifi- 
cant amount of mass. In our current model, we effectively 
keep track of an "inner" disk: the portion containing the 
majority of the mass. This justifies our choice of sur- 
face density profile E oc r _1 , since the radius at which 
this power-law slope is achieved is also the radius that 
encloses most of the mass. We allow for a small amount 
(2%) of material raining in from the core to be carried 
out with the angular momentum. 

The disk's angular momentum is then equal to that of 
the infalling material, in addition to the amount already 
in the disk, minus some portion which has been trans- 
ferred to this outer disk. The disk loses a fraction bj of 
its angular momentum and a small amount of mass on 
the viscous timescale t v = Md/M*, so long as matter is 
still accreting from the core: 



3d — j in Mi 




(22) 



As above, the subscript "in" denotes newly accreted mat- 
ter. The factor (M; n /M*) is roughly unity in the main 
accretion phase, but goes to zero when accretion stops. 
We thus assume the outer disk only applies a torque when 
it gains matter from the inner disk. Without accretion 
the outer disk has no effect, and thus after accretion ends, 
the inner disk is free to expand sclf-similarly at constant 
3d- We consider this a conservative approach, consider- 
ing that we do not treat effects like photoionization that 
might remove material from the inner and outer disk, 
especially in massive stellar clusters. 

We consider bj = 0.5 to be typical; in this case an ac- 
creting disk loses about half its angular momentum each 
viscous time. Since the disk sheds mass at the same rate, 
this allows its radius to remain comparable to the circu- 
larization radius of the infalling material. Although our 
choice of bj is somewhat arbitrary we demonstrate that 
our parametrization makes the disk evolution somewhat 
independent of this value. See §5.6| for discussion. 



3.5. Disk Fragmentation 

Since we have now computed M^, -M*, and J d , our 
model is almost complete. However, as demonstrated 
by both pre vious analytic work (KM06) and numerical 
simulations dKrumholz e t al. 2007b Vorobyov & Basu 
|2006| |Lodato~&; Kice|2005p , our parameter space extends 
deeply into the regime where disk fragmentation is ex- 
pected. We must account for this to model disk evolu- 
tion. It is not our intent to follow the detailed evolution 
of the fragments formed, nor their mass spectrum; we are 
interested primarily in how they help the disk regulate 



In keeping with the approach outlined in § |3.1[ we 
make the important assumption that the disk fragments 
into small objects when Q drops below a criti cal value, 
Qcrit ; which we take to be unity. Other authors ( jGammie 
2001 Rice et aT7||2003 1 have pointed out the importance 



boundary. In particular they find, in simulations with 
imposed cooling, a critical value of T c f2 above which disks 
do not fragment, and below which they do. 

Our fragmentation model reproduces these results (in- 
deed, it is calibrated to the same simulations) and we 
believe that the two views are in fact equivalent. Within 
our model, a disk whose Q is close to unity will be heated 
by accretion at a rate close to the critical cooling rate 
found in these simulations. In the absence of any addi- 
tional heating, the cooling rate must exceed the critical 
value in order for Q to fall below unity, so that frag- 
mentation can commence. In other words, since in our 
model Q is calculated based on the competition between 
cooling and the combination of viscous dissipation and 
irradiation, if Q falls below unity then it is necessarily 
the case that the cooling rate is sufficient to overwhelm 
viscous heating, and therefore to s atisfy a cooling con- 
dition simi l ar to those identified by Gammie (20011 and 
Rice et all (12003 1. 



1'he benefit of our fragmentation model is that it can 
be easily extended into the realistic regime of irradiated 
disks, whereas a model that refers solely to the cooling 
time cannot. 

We note, in support of our model, that we know no 
examples of disks for which Q < 1 that do no t fragment, 
nor those with Q > 1 that do. Moreover, |Rice et al. 



(2003) note that a sufficiently slowly cooling disk reaches 



an equilibrium at a Q value higher than unity; this is 
consistent with a heating rate that drops sharply as Q 
increases, as our accretion model would predict. 

To implement fragmentation within our numerical 
models, we must specify how much mass goes into frag- 
ments each time step when Q < 1. We first define a 
critical density, E dc : 



c„f2 



nGQ, 



crit 



(23) 



a reduction of surface density from E^ to E c would return 
the disk to stability. Because we expect fragmentation to 
happen over a dynamical time, we assume that it depletes 
the disk surface density at the rate: 

E frag = -(E d -E d , c )fi, (24) 

This rate is fast enough to ensure that Q never dips ap- 
preciably below Qcrit- 

For simplicity, we assume that while fragments con- 
tribute to the mass of the disk, they do not enter in 
Toomre's stability parameter Q except insofar as they 
contribute to the binding m ass. (One could consider a 
composite Q: |Rafikov 2001 ) Nor do we follow the mi- 
gration of fragments m the disk. Instead, we allow them 
to accrete onto the central star at the rate 



i'rnjj 



0/-Mf rag fi, 



(25) 



of a disk's thermal physics in setting the fragmentation 



with 4>f = 0.05. The assumption is simply that some 
fraction of the fragments accrete each orbit. Fragments 
form preferentially at large distances from the star, and 
thus only a small amount of the fragment mass will make 
it into the central star each orbit. Changing this param- 
eter by an order of magnitude only marginally alters the 
disk evolution. 

We also make the important assumption that disks will 
always fragment to maintain stability, and allow accre- 
tion to proceed. While this is likely a good assumption 
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based on the existence of massive stars that appear to 
have formed via disk accretion, the persistence of rapid 
accretion during fragmentation has not been satisfacto- 



rily demonstrated in numerical simulations. See ^7.1 



3.6. Binary Formation 

A majority of stars, especially massive stars, are found 
in binary and multiple systems. Though we present a 
very simplified scenario for star formation, we do account 
for the possibility that a single secondary star will form 
if Md > M*, that is, if the disk grows unphysically l arge 



with respect to the central star. (As we discussed in § 3.4 
this may well be conservative - in the sense that secon- 
daries may form at even lower values of /i, or at earlier 
times through c ore fragmentation as described in |Bon-| 
nell et al. 2004 ) When [i > 0.5, we remove the excess 
mass and store it (and the associated angular momen- 
tum) in a binary star. Because this tends to happen 
before the disks have become very extended, we assume 
the binary separation will be small; we therefore ignore 
the binary as a source of angular momentum for the disk. 
As with fragments, we assume the disk is affected by the 
secondary star only through the increased binding mass. 
We make no attempt to account for its contribution to 
the total luminosity. 

3.7. Summary of Model 

We summarize our model via the flowchart shown in 
Figure [3] which illustrates a simplified version of the 
code's decision tree. At a given time t we know the 
current disk and star mass, and the current angular mo- 
mentum and mass infall rates as prescribed in §2.1 and 



{2.2 We can calculate Rd and directly, and find the 
appropriate stellar luminosity based on its evolution, cur- 
rent mass, and accretion rate. Using these variables we 
self-consistently solve for the appropriate temperature, 
Q, and disk accretion rate as described in {3.3 With 
this information in hand, we determine whether the disk 
is stable, locally fragmenting, or forming a binary. If 
the disk is stable, we proceed to the next iteration. If 
Q < 1 , then the disk puts mass into fragments according 
to equation [23] . If /i > 0.5 we consider binary forma- 
tion to have occurred, and the net angular momentum 
and disk mass over the critical threshold is placed into 
a binary (see { 3.6 1. We stop simulations after 2 Myr 
for two reasons: first, the most massive stars in our pa- 
rameter space are significantly evolved and so our stellar 
evolution models arc no longer sufficient; and second, 
because many other effects begin to dominate the disks 
appearance at late stages due t o gas-dust interaction and 
photo-evaporation (Keto 2007). 



4. EXPECTED TRENDS 

Before examining the numerical evolution, it is useful 
to make a couple analytical predictions for comparison. 

First, can we constrain where disks ought to wander in 
the plane of Q and /i? This turns out to depend critically 
on the dimensionless system accretion rate 



so. = 



MhJL 



M* d Q(R cixc ) G 2 M* d 



(26) 



which is the ratio of the mass accreted per radian of disk 
rotation (at the circularization radius i? c irc) to the total 



system mass M*<j = M* + M d - Since the active inner 
disk has a radius comparable to i? c irc, this controls how 
rapidly the disk gains mass via infall. 

The importance of 5ft; n is apparent in the equation gov- 
erning the evolution of the disk mass ratio \x: 



A = M in fl_ _ 

= gCgclrc) ( 1 



M, 



1 Ki, 



MM' 



(27) 



Since we consider M* / (Mdty to be a function of /i and 
Q, we must know the dis k temperature to solve for fj.(t). 
Regardless, equation (27 1 shows that larger values of 



tend to cause the disk mass to increase as a fraction of 
the total mass. We may therefore view 5Rj n and Q as the 
two parameters that define disk evolution - of which 5R; n 
is imposed externally and Q is determined locally. 

Moreover, 3?i n takes characteristic values in broad 
classes of accretion flows, such as the turbulent core mod- 
els we employ. Suppose the rotational speed in the pre- 
collapse region is a fraction f K of the Kepler speed, so 
that j in = f K rv K (r) = /^[GrAf^r)] 1 / 2 , and suppose 
also that the mass accretion rate is a fraction e/ acc of 
the characteristic rate VK(r) 3 /G. Then, 



so. — 

J tin 



f K /ace 



(28) 



(In this expression, negative three powers of e appear 
because the binding mass is e times smaller for the disk 
than for the core; one of these is compensated by the 
reduction of the accretion rate.) 
In 



2.1 



_ we adopted the |McKee fc Tan| ( |2003| model for 
massive star formation due to core collapse of a singular, 
turbulent, polytropic sphere in initial equilibrium. Their 
equations (28), (35), and (36) imply 



/ acc = 0.84(1 - 0.30fc p ) 



3 hp 
l + H 



1/2 



(29) 



within 2%, where 1 + Hq ~ 2 represents th e support due 
to static magnetic fields ( |Li fc S hu 1996). (Note , their 
equation [28] is a fit made p y McKcc & Tan 2002 to the 
results of McLaughlin fc Pudritz||1997| ) 

KM06 predicted the turbulent angular momentum of 
these cores; our parameter fx equals {Oj^j) 1 ' 2 in their 
paper. Their equations (25), (26), and (29) imply 



!k = 



0.49 (1 - k p /2f 

,1/2 



(fc„-l)V2 



(30) 



with excursions upward by about 50% and downward by 
about a factor of three expected around this value; here 
4>b ~ 2.8 represents the magnetic enhancement of the 
turbulent pressure. All together, we predict 

1.26 



0.10 

q.3/2 



■0.02 



l + #o 



1/2M- 



-(1 -0.30fc p ) 



where the evaluation uses 1 + Hn 



(31) 



and 



fan 



1.5. 
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Core Model 

M in (t),3 in (t) 



stop ) 

..♦ yes 




fragments 

2d |, Afftag T 



start 



Initial Conditions 

M* , M d0 , J d0 




global disk properties 

M d , J d , M* 



Thermal Physics 

► viscous heating 

► irradiation 
^ radiative cooling 




Disk Braking Model 



binary 

M d l,J d l,M bin T 



Fig. 3. — A simplified schematic of the decision tree in the code. The primitive variables, M<j, M*, and Jj, together with the core 
model, Mi n (t) and Ji n (t), allow for the determination of all disk parameters at each time step. Note that c B ,Q, and M are solved for 
simultaneously. Once t he se lf-consistent state is found, the values of Q and fi determine whether either the binary or fragmentation regime 
has been reached. See j]3.7|for a description of the elements in detail. 



Importantly, 5ii n is a function of (1 + Hq), cf>B, e, and 
k p , but not the core mass. We therefore expect similar 
values of 5fj n to describe all of present-day massive star 
formation, at least insofar as these other parameters take 
similar values. Suppose, for instance, that the formation 
of 10M Q and 100M Q stars were both described by the 
same 3?i n . According to equation (27), the difference in /x 
between these two systems would be controlled entirely 
by the thermal effects that cause them to take different 
values of Q. 

A few additional expectations regarding Q itself can 
be gle aned from the analytical work of Matzn er fc Levin] 
(2005) and KM06: 



The Toomre parameter remains higher than unity 
for low-mass stars (< IM@) in low-column cores 
(E Cj o <C 1), but falls to unity during accretion for 
massive stars and for low-mass stars in high-column 
cores; 

A given disk's Q drops during accretion, reach- 
ing unity when the disk extends to radii beyond 
~ 150 AU (in the massive-star case), or to peri- 
ods larger than ^460 yr (in the case of an optically 
thick disk accreting from a low-mass, thermal core) . 



- At the very high accretion rates characteristic of 
the formation of very massive stars (> 1.7 x 
lO~ 3 Af yr" 1 ), disk accretion is strongly destabi- 
lized by a sharp, temperature-dependent drop in 
the Rosseland opacity of dust. 

(For more detailed conclusions, see the discussion sur- 
rounding equation [35] in (KM06).) With the help of 
equation (27) we also deduce that more massive stars 



will have generally higher disk mass fractions, because: 
(1) they are described (in our model) by the same value 
of 3f?i n ; (2) more massive stars achieve lower values of Q; 
and (3) in our model, lower Q leads to lower values of 
M*/(M(jf2), so long as Q > 1.3. The conclusion that 
higher-mass stars have relatively more massive disks fol- 
lows from these three points by virtue of equation |27] . 

More generally, any effect which causes M* / '(M^Q) to 
drop (without affecting 3?i n ) will tend to increase /it, and 
vice versa; this conclusion is not limited to our adopted 
disk accretion model. 

Within our model, JW*/(Mdf2) increases with Q unless 
1 < Q < 1.3, in which case the dependence is reversed. 
Disks ought therefore to traverse from high Q and low 
fi, to low Q and high /i, until Q = 1.3; but for 1 < 
Q < 1.3, fj, and Q should decline together. In physical 
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terms, this reversal represents a flushing of the disk due 
to the strong angular momentum transport induced by 
the short wavelength gravitational instability. 

We now turn to our suite of numerical models to test 
these expectations. 

5. RESULTS 

We begin by examining the evolution of disks through 
their accretion history for a range of stellar masses, de- 
termining when and if they are globally or locally unsta- 
ble, and the dominant mechanism for matter and angu- 
lar momentum transport through disk lifetimes. Next, 
we explore how these results are influenced by varying 
the other main physical parameters: T c , S c , omri, and 
by varying the angular momentum prescription. For this 
purpose we first define a fiducial sequence of models in 
§ |5.1[ and then expand our discussion to the wider param- 
eter space encompassed by the aforementioned variables. 
Figures |4] - [6] show results from our fiducial model, and 
Figures [8] and [9] explore the effects of our environmental 
variables. 

Because our prescription for disk accretion and frag- 
mentation is necessarily approximate, any specific pre- 
dictions are unlikely to be accurate. We concentrate in- 
stead on drawing useful observational predictions from 
our models' evolutionary trends. 

5.1. The fiducial model 

Our fiducial model explores a range of masses with a 
standard set of parameters, which we list in Table [T] For 
our exploration of the stellar mass parameter space, we 
allow E c to vary as E c = 10" 1 - 84 (M C /M Q ) - 75 (with an 
enforced minimum at 0.03 g cm~2 so that S c varies from 
0.03 — lg cm -2 across the mass range 0.5 — 120M Q ). This 
relationship ensures that for our fiducial model, each sys- 
tem is forming at a S c that is characteristic of observed 
cores. Enforcing this S c — M c correspondence specifies 
the core r adiu s. We explore the effects of S c indepen- 
dently in £5.3 All "low mass" runs that are shown e.g. 
1M Q , have X c ,iow = 0.03 g cm" 2 , and "high mass", e.g. 
15M Q , have S Cihigh = 0.5 g cm" 2 . All systems start 
out with an initial stellar mass of 0.10Af Q , disk mass of 
0.01M Q and jd — 10 19 . Varying these parameters over 
an order of magnitude effects the initial evolution for a 
few thousand years, but runs converge quickly. One can 
find pathological initial conditions, particularly for small 
mass values. We believe this is due to the lack of sen- 
sitivity of a one zone model. The initial disk radius is 
calculated self-consistcntly from the amount of mass col- 
lapsed into the system at the first time step, the initial 
Jd is typically smaller by a factor of a few than ji n . 

As illustrated by the evolutionary tracks of accreting 
stars in the Q — fi plane in Figure [4] our model agrees 
with the general trends of previous work and with the 
expectations described in §|4] in that low mass systems 
are stable and have low values of /i, while more massive 
syste ms undergo a period of strong gravitational insta- 



bility (Krumholz ct al. 



2007b', KM06). Here we see that 



as we go to higher stellar masses, disks spend more of 
their time at high /i and undergoing disk fragmentation. 
For stars of < 1M Q , Q stays above unity, and the disks 
remain Toomre stable, although still subject to gravita- 
tional instability due to their non-negligable disk masses 
(see Figure [6]) . (Note that due to our abrupt shift in the 



Parameter 


Fiducial 


Range 


3 


0.5 


- 1.0 


^c,low 


0.03 g cm" 2 


0.03 - 1 g /cm 2 


^c,high 


0.5 g cm -2 


0.03 - 1 g /cm 2 


T c 


20 K 


10 - 50K 


Q MRI 


0.01 


0.001 - 0.1 


4>f 


0.05 


- 0.5 


e 


0.5 


N/A 



TABLE 1 

Fiducial parameters for disk models for low and high mass 
stars, and the accompanying ranges explored. 



disk irradiation model, there is a small discontinuity in 
the temperature calculation at the end of accretion which 
can cause unphysical fragmentation even at low masses, 
and a jump in Q at all masses.) The expectation that Q 
and \i evolve in opposite directions until Q < 1.3 is also 
roughly borne out. However, note that for the 15M Q 
star-disk system (right plot), the accretion rate is great 
enough that there is a build up of mass in the disk once 
Q reaches unity, and the local instability s atur ates. This 
saturation leads to binary formation (see {5.7). 

Figure [5] shows the evolution of Q through the accre- 
tion history of a range of masses. We see that disks be- 
come increasingly susceptible to fragmentation with in- 
creasing mass. Disks born from cores that are smaller 
than about 2Mq remain stable against fragmentation 
throughout their evolution, although we expe ct moder- 
ate s piral s tructure (as is seen in the models of Lodato & 
20041. Recall that with e = 0.5, a 2M Q core makes 



Rice 



a IMq star-disk system. Figure [6] illustrates the corre- 
sponding evolution of /i throughout the accretion history 
for the same set of masses. As described in Qthe typical 
disk mass increases with stellar mass. At high masses, 
binary formation occurs during the peak of accretion just 
before 10 5 years, and for stars > 100M Q , there is an early 
epoch of binary formation at roughly 10 4 years. 

We also see that for higher mass cores there are three 
relatively distinct phases through which disks evolve: 

- Type I: Young, < 10 4 yr systems, whose disks 
are described by small mass fractions and rela- 
tively high Q. These would appear as early Class 
sources, deeply embedded in their natal clouds. 

- Type II: Systems between 10 4 — 10 5 yrs in age, 
whose disks are subject to spiral structure, and in 
high mass systems, fragmentation. Disk mass frac- 
tions are ~ 30% — 40%, substantially higher than 
in Type I systems. These disks would appear in 
Class 0-1 sources. 

- Type III: Systems older than 10 5 yrs, which have 
stopped accreting from the core, and consequently 
acquire low disk mass fractions as the disks drain 
away. These are the disks that are most like those 
observed in regions of LMSF as Class I objects. 

These three stages serve as a useful prediction for future 
observations; see <|6]for more details. 

5.2. Influence of vector angular momentum 

The accretion disk's radius plays a critical role in deter- 
mining whether or not the disk fragments. Consequently, 
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Fig. 4. — Evolutionary tracks in the Q, /i plane of a lMQ(left) and 15 Mq (right) final star-disk system overlayed on the contours of our 
accretion model (contour spacing is identical to Figure |2t. The white arrows superposed on the tracks show the direction of evolution in 
time. The low mass star remains stable against fragmentation throughout its history, while the more massive star undergoes fragmentation 
and more violent variation in disk mass. The jump a the end of accretion in the 15 Mq system is due to the switch in the irradiation 
calculation. 




years 

Fig. 5. — Contours of Q over the accretion history of a range 
of masses for the fiducial sequence. Masses listed on the y-axis 
are for the total star-plus-disk system final mass — because the 
models halt at 2 Myr, some mass does remain in the disk. Contours 
are spaced by 0.3 dex. At low final stellar masses, disks remain 
stable against the local instability throughout accretion. At higher 
masses, all undergo a phase of fragmentation. One can see three 
distinct phases in the evolution as described in j|5] Disks start 
out stable, subsequently develop spiral structure as the disk mass 
grows and become unstable to fragmentation for sufficiently high 
masses. As accretion from the core halts, they drain onto the star 
and once again become stable. 



we expect our results to depend somewhat on effects that 
change the disk's specific angular momentum. Because 
we track the vector angular momentum of the inner disk, 
and because our turbulent velocity field is three dimen- 
sional, we account for a possible misalignment between 
the disk's angular momentum axis, J, and that of the 
infalling angular momentum, j; n . The wandering and 




Fig. 6. — Contours showing the evolution of fi = M for 

the fiducial sequence. Each contour shows an increase of 0.05 in 
fx. Again one can see the division into three regimes: low mass 
disks at early times, higher mass, unstable disks that may form 
binaries during peak accretion times, and low mass disks that drain 
following the cessation of infall. Systems destined to accrete up to 
~ 70A/q or more experience two epochs of binary formation in our 
model. In these systems the accretion from the core exceeds the 
maximum disk accretion rate very early, causing the disk mass to 
build up quickly. 

partial cancellation that result provide a more realistic 
scenario than given by the KM06 analytic approxima- 
tions, in which vector cancellation is accounted for only 
in an average sense. In practice, however, the disk and 
infall remain aligned rather well (J • j; n ~ 0.8), so mis- 
alignment plays only a minor role in limiting the disk 
size. This is illustrated by Figure [7j in which we com- 
pare the disk radius in two numerical realizations of the 
turbulent velocity field, against one in which j m has a 
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-fiducial model 
-alternate realization 
-KM06 
eircularization radius 




Fig. 7. — Comparison of disk radius over the evolution of a 20 
Mq star-disk system in four cases: the KM06 analytic calculation, 
the eircularization radius of the currently accreting material, and 
two realizations of the numerical model. The analytic case over- 
estimates the expected radius at early times because it does not 
allow for cancellation of vector angular momentum. Similarly the 
eircularization radius is an overestimate because the disk has no 
"memory" of differently oriented j. At later times, the eircular- 
ization radius approaches the standard radius calculation for that 
realization (thick black line) demonstrating the concentration of 
turbulent power at large scales. 

fixed direction and obeys the KM06 formulae. We also 
plot the infall eircularization radius i? c irc (of one numer- 
ical realization) for comparison. In general we find that 
the analytic prescription slightly over-predicts the disk 
radius at early times; this is partly due to "cosmic" vari- 
ance in the numerical realization, and partly due to disk- 
infall misalignment. 

5.3. Varying E c 

We explore the effect of individual parameters by con- 
sidering one or two systems along our fiducial sequence, 
and varying parameters one by one relative to their fidu- 
cial values. First, we vary S c over 0.03-1 gem -3 , span- 
ning the r ange from isolated t o intensely clustered star 



formation (Plume et al. 1997). Column density affects 



star formation in two primary ways: it influences the core 
radius (by determining the confining pressure) and the 
accretion rate during collapse (again, by setting the outer 
pressure and thus the velocity dispersion). However, 
these two effects counteract one another: smaller values 
of S c correspond to larger cores and larger, thus more un- 
stable disks (Rd oc E c ' ), but smaller E c also leads to 

lower accretion rates and thus stabler disks (M cx S^ 4 ) . 
The thermal balance of the disk midplane is affected by 
these trends. An analysis by KM06 (see their equation 
[35]) shows that higher S c inhibits fragmentation if the 
disk temperature is dominated by viscous heating (which 
is proportional to the accretion rate), but enhances frag- 
mentation if irradiation dominates (when the increase in 
accretion generated heating is insignificant), and that the 
two effects are comparable along our fiducial model se- 
quence. We therefore expect fragmentation to be quite 
insensitive to E c , for massive star formation along our 
fiducial sequence. This is precisely what we find in our 
models: disks born from lower-S c cores, in lower pres- 
sure environments, evolve in essentially the same way, 
but more slowly. 
In contrast, disks around low-mass stars - those with 



final masses comparable to the ther mal Jeans mass - 
are st able at low E c (as predicted by Matzner & Levin 
|2005 1, and because irradiation dominates at larger radii, 
higher S c tend to enhance fragmentation there. Figure 
M illustrates the evolution of Q for a IMq accreting star 
for a range of column densities. 

5.4. Varying T c 

Observations of infrared dark clouds, and sub-mm 
core detections find typ ical temperatures from 10 — 50K 
(Johnstone et al. 20011. In our models T c determines 
the amount of thermal, and therefore turbulent, support: 
higher temperatures require less turbulent support in the 
core. Temp erature also sets the th ermal Jeans mass M t ^ 
within the McKee & Tan ( 2003 1 two component core 
model. Accretion from this thermal region leads to more 
stable disks; therefore, higher core temperatures increase 
the mass of a star which can accrete stably. 

Figure [8] shows the evolution of Q during the accretion 
of a system with final mass 1M© over a range of temper- 
atures (all other parameters take their fiducial values). 
The difference in evolution is negligible for high mass 
stars, as these accrete from supersonic cores. 

5.5. Varying q;mri 

This work is not an exploration of the detailed behav- 
ior of the MRI; we include it as the standard mechanism 
for accretion in the absence of gravitational instabilities, 
which in most sc enarios (asi de perhaps from low mass 
stars whose disks Shu et al. 2007 have argued may be 
strongly sub-Keplerian) overpower the MRI. However, 
the strength of the MRI does influence the transition to 
gravitationally dominated accretion in the Q — fj, plane 
as shown in Figure [2] The strength of the MRI also 
influences the maximum disk mass obtained before grav- 
itational instabilities set in: higher values of «mri re- 
duce the influence of gravitational instabilities by in- 
suring that the disk drains more quickly, whereas lower 
values expedite the transition to gravitational instability 
driven accretion. Figure [9] shows the influence of aMRi 
on jtx; the influence on Q is less dramatic: the descent 
of Q towards unity is marginally delayed for the strong 
MRI case. 

5.6. Varying bj 

Our most uncertain variable is the braking index bj, 
which determines the rate of angular momentum ex- 
change with an outer disk. However, disk evolution turns 
out to be rather insensitive to this parameter. The pri- 
mary reason for this is the concentration of power in the 
turbulent velocity field on the largest scales: ji n is always 
large compared to the disk-average j. This reduces the 



importance of the loss term in equation (22 1. As a result, 
although the period in which the disk is fragmenting is 
reduced in the high bj case, it is merely postponed by 
~ 10 4 yrs. The braking index does have moderate in- 
fluence on the disk mass during the peak of accretion, 
and thus on the formation of binaries. Figure 9] shows 
the evolution of \x for a system accreting towards 15 Mq 
from a 30M Q core. Here one can see the influence of 
bj on binary formation. Low values of bj correspond- 
ing to higher net angular momentum produce binaries at 
lower masses by allowing the disk mass to grow larger. 



13 



30 



r- 



10 L 



10° 



10' 



10" 



10" 




0.5 
0.4 

or 

o 

0.3 | 

0.2 

0.1 



years 



10 10 

years 



Fig. 8. — Contours of Q showing the effect of initial core temperature T c (left) and S c (right) on the evolution of a 1 Mq final star-disk 
system. Contour spacing is 0.1 dex (except the lowest two contours which are spaced by .05 dex). Increasing S c tends to marginally 
destabilize the disk, while higher temperatures stabilize the disk. We exclude temperatures too high for the 2Mq core to collapse given its 
initial density, i.e., those above 40 K. 



Notably, even for the 15M Q final mass star shown here, 
the smallest mass for which binaries form in the fiducial 
model, the change in disk mass is only ~ 10%. 

5.7. The formation of binaries 

Within the context of our model for disk fission into 
a binary system, (described in 



3.6 1, the formation of a 



companion is strongly dependent on the infalling angular 
momentum distribution, and on the turbulent velocity 
profile of the particular core. In our fiducial model, bi- 
nary formation occurs for cores above 3OM . For cores 
> 140Mq, there are two epochs of binary formation, the 
first one at roughly 10 4 years. This mass boundary is 
quite sensitive to our conservative threshold for disk fis- 
sion, yU = 0.5: binaries may well form at lower values of /i, 
and thus at lower masses (see Figure [6]) . The mass of the 
binaries that form increases with initial core mass. This 
increase simply indicates that the mass ratio exceeds the 
critical value for more time, as we do not include a mech- 
anism for accretion onto the binary. As such, we do not 
predict values for the binary mass ratio q, but simply 
indicate the regimes in which binary formation seems 
likely. The 30M Q core cut-off is fairly robust to varia- 
tions in S C ,T C and bj over the ranges discussed above 
for our fiducial turbulent field. Cosmic variance in the 
field from one realization to another has a much larger 
effect on binary formation than any of our other model 
parameters (aside from /i C rit)- 

In our fiducial model disk fission only occurs when the 
gravitational instability has saturated and Q ~ 1. This 
means that the disk is draining at the maximum rate 
given its mass. If matter is falling in from the core more 
rapidly than this rate, the disk mass will increase: if 
the accretion rate from the core exceeds the maximum 
rate at Q = 1 and fx = 0.50, disk fission occurs. In our 
fiducial model, this corresponds to an accretion rate on 
to the disk: M in /M d tt = 10~ 2 - 36 . The early epoch of 
binary formation at very high masses is a consequence of 



this limit: since the accretion rates begin to exceed the 
critical rate sooner, the disk's mass increases earlier in 
its evolution. This critical value is in agreement with the 
prediction of KM06 that disks are sharply destabilized 
when accreting at rates higher than 1.7 x lO~ 3 A/ yr _1 . 

The time at which binaries form is also very depen- 
dent on the angular momentum profile. In the fiducial 
model, the lowest mass binaries form during the peak of 
accretion, at ~ 10 5 years, but as the final system mass 
increases, binary formation pushes to earlier times ~ 10 4 
years. In certain runs we find earlier binary formation at 
smaller masses (< 10 4 years) when there is a peak in the 
infalling angular momentum profile which rapidly sends 
Q towards unity. The presence of binaries in much of 
our parameter space illustrates that heavy circumbinary 
disks may be critical to binary evolution. 

Observations suggest that a range of binary systems 
exist as a result of variations in angular momentum as 
evidenced by the presence or lack of disks around each 
component. Submillimeter observations of lower mass 
objects in Taurus have revealed evidence for a binary 
with circumstellar and circumbinary disks (Osorio et al 



2003 1 , where the binar ies are close enough t o cause disk 
truncation (~ 45 AU). Anglada et al. (20041 have found 



another Class 0/1 binary system in 1NGC 1333 in which 
only the primary has a disk: the diversity of systems is 
likely due to the variat ions in angular momentu m of the 
infalling material. As |Bate fc Bonncll (19971 suggest, 
binaries forming from low angular momentum material 
will likely not form their own disks, while those with 
higher angular momentum may. It seems plausible that 
the absence or presence of secondary disks is indicative 
of the formation process of the system. 

As illustrated by these observations, the dependence 
we find on core angular momentum is a sensible out- 
come: one expects the chance rotation to have a stronger 
effect on multiplicity than other parameters like temper- 
ature and density, which set the minimum fragmenta- 
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Fig. 9. — Contours in fi illustrating the influence of varying omri 
(top) and the braking index, bj (bottom) for a star-disk system of 
final mass 15Mq, the lowest mass at which a binary forms in our 
fiducial model. Contours of fi are spaced by 0.05. The upper plot 
shows the effect of varying cimri from 10 5 — 10 -1 ' 5 . While the 
change has little effect on the evolution of Q, the disk fraction fi 
decreases with increasing «mri- As a result, the mass at which 
binary formation begins is pushed to higher masses. The lower 
plot shows the effect of varying the braking index bj . An increase 
of bj lowers the disk angular momentum, reducing the disk mass 
and inhibiting binary formation. Note that the variation in disk 
mass is only ~ 10%. 

tion mass. We emphasize that we are only exploring 
one possible path for binary formation, and predict that 
disk fragmentation is an important, if not the dominant 
mechanism at high masses and colu mn densitie s. This 



is especially true since, as argued by Krumholz (20061, 
who has shown that for massive stars, once the central 
core has turned on, the Jeans mass rapidly increases due 
to the stellar luminosity, significantly reducing the pos- 
sibility for Jeans-instability induced core fragmentation. 

6. OBSERVABLE PREDICTIONS 

Our models make strong predictions for the masses 
and morphologies of disks during the embedded, accret- 
ing phase, and these will be directly testable with future 
observations. Detailed calculations based on radiation- 
hydrodynamic simulations of massive protostellar disks 
indicate that disks with fi of a few tenths around stars 
with masses > 8 M Q , corresponding to embedded sources 
with bolometric luminosities > 10 4 Lq, should produce 
levels of molecular line emission that are detectable and 
resolvable with ALMA in the sub-millimeter out to dis- 
tances of a few kpc, and with the EVLA at centimeter 
wavelengths at distances up to ~ 0.5 kpc (Krumholz 



par- 
ticularly efficient at observing protostellar disks, since 
ALMA's large collecting area will enable it to map a mas- 
sive disk at high resolution in a matter of hours. Dust 
continuum emission at similar wavelengths should be de- 
tectable at considerably larger distances, although the 
lack of kinematic information associated with such ob- 
servations makes interpretation more complex. Regard- 
less of whether dust or lines are used, observations using 
ALMA should be able to observe a sample of hundreds 
of protostellar disks around embedded, still-accreting 
sources, with masses up to several tens of M®. 

The main observational prediction of our model is the 
existence of type II disks - those with /i of a few tenths or 
greater and Q«l - and the mass and time-dependence 
of the type II phase. Examining Figures [5] and [6j we 
see that our model predicts that protostellar cores with 
masses < 2 M Q should experience only a very short type 
II disk phase, or none at all. In contrast, cores with 
larger masses have type II disks for a fraction of their 
total evolutionary time that gets larger and larger as the 
core mass rises, reaching the point where type II disks 
are present during essentially the entire class 0, accreting 
phase for cores > 100M Q in mass. 

Type II disks have several distinct features that should 
allow observations to distinguish them from type I or 
type III disks, and from older disks like those around T 
Tauri and Herbig AE stars. First, since type II disks are 
subject to strong gravitational instability, they should 
have strong spiral arms, with most of the power in the 
m = 1 or m = 2 modes. This is perhaps the easiest 
feature to pick out in surveys, since it simply requires 
observing the disk morphology and can therefore be mea- 
sured using continuum rather than lines. 

Second, because their self-gravity is significant, type II 
disks will deviate from Kcplerian rotation due to non- 
axisymmetric motions, and will also be super-Keplerian 
in their outer parts compared to their inner ones. The 
latter effect arises because, when the disk mass is com- 
parable to the stellar mass, the enclosed mas s rises as 
one m o ves ou tward in the disk. Recent work by |Torrelles] 



et al. (20O7J provides a possible example of this phe- 
nomenon^ The source HW2 in Cepheus A is predicted 
to have a central mass of order 15M , and a disk ra- 
dius of 300 AU, with a temperature s lightly under 200 K 

2007). High resolution 



(|Patel et al.|2005| |Torrelles et al. 



VLA obs ervations now show eviden ce of non-Keplcrian 



rotation ( |Jimenez-Serra et al.|2007 l, consistent with our 
predictions for type 11 disks. 

Third, a type II disk is massive enough for the star- 
disk system center of mass to be significantly outside of 
the stellar surface if the disk possesses significant non- 
axisymmetry. As a result, the star will orbit the center 
of mass of the system, and this will produce a velocity 
offset of a few km s _1 between the stellar velocity and 
the zero velocity of the inner, Keplerian parts of the disk 



Bertin fc Lodato|1999| |Rice et al.|2003|lKrumholz et al 
2007ap . This should be detectable if the stellar velocity 
can be measured, which may be possible using Doppler 
shifts of radio recombination lines for stars producing 
hypercompact HII regions, or using proper motions for 



stars with large non-thermal r adio emission (Bower et al 
2003 1 . In fact recent work by Torrelles et alT j2007| has 



observed said offset. As suggested by |Lodato fe Bertin| 
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(2001 |2003 1, one could also look for the effect in the 
unresolved radio emission from FU Orionis objects. 

A final point concerns the limited range of the disk-to- 
system mass ratio in our simulations, with 0.2 < // < 0.5 
during most of embedded accretion (our type II disks) . 
The upper envelope of fi depends in part on our binary 
fragmentation threshold /i = 0.5. However, in the ab- 
sence of disk fission, disks in our fiducial model never 
grow larger than [i = 0.55. The fact that most accre- 
tion occurs with /i ~ 0.3 provides strong evidence that 
accretion disks do not become very massive compared t o 



the central point mass (as argued by Adams ct al. 1989) 
Current observations such as those of[Cesaroni|(|2005p de 



scribe massive tori with sub-Keplerian rotation and com- 
parable infall and rotation velocities. These structures 
are distinct from the disks that we model: our finding 
that disks hover around \i = 0.3 suggests that higher res- 
olution observations may reveal the Keplerian structures 
within the tori. The underlying physical reason for this 
is that it is not possible to support a mass comparable to 
the central star in a rotationally supported disk for long 
periods of time; gravitational instabilities will destabilize 
such a disk on orbital timescales, causing it to lose mass 
either through rapid accretion or fragmentation. 

7. CONCLUSIONS 

We have constructed a simple, semi-analytic one-zone 
model to map out the parameter space of disks in Q — /j, 
space across a range of stellar masses throughout the 
Class and Class I stage, pushing into the Class II phase. 
We include angular momentum transport driven by two 
different mechanisms: gravitational instability and MRI 
transport modeled by a constant a. Our model for an- 
gular momentum infall is unique in that we keep track of 
an inner and outer disk, and infall direction so that can- 
cellation may occur as the infall vector rotates. We allow 
for heating by the central star, viscous dissipation and 
a background heat bath from the cloud accounting for 
both the optically thin and optically thick limit within 
the disk and accreting envelope. By requiring that the 
disk maintain mechanical and thermal equilibrium, we 
determine the midplanc temperature at each time step, 
and thus Q in the disk. 

7.1. Caveats 

In interpreting the results of our calculations, it is im- 
portant to keep several caveats in mind. Our model 
for fragmentation, though rooted in simulations, includes 
one important assumption: no matter how violently un- 
stable a disk becomes, it can always fragment, return to 
a marginally stable state, and continue accreting. While 
the existence of stars well into the mass regime of frag- 
mentation makes this outcome seem likely, it has yet to 
be demonstrated in simulations. Equally untested is the 
hypothesis that when fragmentation is strong enough , 
i.e., when M in > c 3 JG so that Q < 1 ( |Gammie||200l| , 
accretion onto the central star will be choked off. KM 06 
have argued that accretion is sharply destabilized when 
its rate exceeds 1.7 x 10~ 3 Af Q yr _1 , due to a drop in the 
Rosseland opacity, and that this may be related to the 
stellar upper mass limit. 

In order to explore a wide parameter space, we do 
not carry out detailed hydrodynamic calculations to de- 
termine the onset of instability, but instead use results 



from previous simulations, and develop analytic formulae 
that describe behavior intermediate between the regimes 
which they explore. Although our approach is very ap- 
proximate, it can be made increasingly more realistic as 
additional numerical simulations become available. Due 
to our one-zone prescription, we cannot resolve spiral 
structure or measure the degree of non-Keplerian motion. 
In addition, we do not follow the evolution of fragments, 
nor their interaction with the disk. Although we allow for 
the formation of binaries, we do not follow their evolution 
and accretion, which limits our ability to make predic- 
tions about mass ratios and angular momentum transfer 
between the disk and the companion. Our model for 
angular momentum infall is responsible for the largest 
uncertainty in our conclusions because different realiza- 
tions of the turbulent velocity field can alter the disk size 
at a given epoch by a factor of a few. Nevertheless, these 
variations are well within the analytic expectations for 
range of angular momenta in cores (KM06). Moreover, 
our approach aims only to predict characteristics of the 
outer accretion disk, and lacks the resolution to track the 
radial profiles of the disk's properties. 

Lastly, recall tha t our models rely on the fundamental 
assumption (§ 3.1 1 that a disk's behavior can be sep- 



arated into dynamical and thermal properties, and in 
particular that its dynamics are governed primarily by 
its mass fraction /i and Toomre parameter Q. 

With these caveats in mind, we summarize our results 
for two different regimes: < 2M Q and > 2Mq. 

7.2. The Low Mass Regime 

Our fiducial models predict that low mass stars will 
have higher values of /x than typically assumed during 
early phases of formation. However, they should remain 
stable against fragmentation throughout their evolution, 
dominated by MRI, long wavelength gravitational in- 
stability, and once again MRI through their evolution 
through the three types of disks discussed in [j5] Dur- 
ing the main accretion phase, disks will have masses of 
order 30% of the system mass. Typical outer radii are 
of order 50 AU, with outer temperatures of 40 K dur- 
ing the main accretion phase, dropping to ~ 10 K at 2 
Myr. The surface density is 10 — 20 g cm~ 2 during the 
main accretion phase, dropping off rapidly at late times 
causing the disk to become optically thin to its own ra- 
diation. As accretion shuts down, and disks grow due to 
conservation of angular momentum, two key effects must 
be considered: truncation and heating by other stars. At 
distances of 1000 AU, very tenuous disks are prone to 
truncation by passing stars particularly in denser clusters 
where average stellar densities are as high as 10 5 stars 
pc -3 ( Hillenbrand fc Hartmann||l998[ ). Similarly, as the 
disk edge extends towards other, potentially more lumi- 
nous stars, the actual flux received will increase, heating 
the d isk above the ~ 10K temperature that we routinely 
find ( | Adams et aT]|2006[ ). 

For core column densities more typical of high-mass 
star forming regions, local instabilities do set in, despite 
the stabilizing influence higher temperatures associated 
with these regions (neglecting the effects of nearby stars) . 
This implies that environment may be important in un- 
derstanding disk evolution 



In contrast to our previous work (Kratter & Matzner 
2006 Matzner & Levin 20051, we find fragmentation at 



1G 



smaller radii. This is primarily due to our modified model 
for aGi, which predicts lower accretion rates and con- 
sequently more fragmentation then previously assumed. 
We note that our results for low-mass systems (final mass 
~ IMq) are rather sensitive to details of the model, such 
as the value of «mri and the way it is combined with aci ■ 



7.3. The High Mass Regime 
For more massive stars, we find high values of /i 



0.35 



and an extended period of local fragmentation as the ac- 
cretion rates peak. Temperatures at the disk outer edge 
at ~ 200 AU approach 100K for systems > 15M Q dur- 
ing accretion, surface densities hover around 50 g cm~ 2 
during the main accretion phase, although by 2 Myr, the 
disks become optically thin in the FIR, as expected. Bi- 
nary formation occurs regularly for cores of order 30M Q 
and higher, though as discussed in §3.6| this is strongly 
dependent on the cosmic variance of the angular mo- 
mentum: cores as small as 20M Q form binaries in our 
model when there is excess angular momentum infall. 
Although fragments accrete with the disk according to 
equation [25] , more massive stars maintain a small mass 
in fragments (10 _1 — 1O _2 M0) in the disk when we end 
our simulations, suggesting that fragments may persist 



to form low mass companions, as predicted in KM06 and 
suggested by the simulations of Krumholz et al. ( 2007b ) . 

Unlike their low mass counterparts, the conclusions we 
draw for massive stars are minimally effected by the en- 
vironmental variables in our model. For the entire range 
of temperatures, densities, and nearly all angular mo- 
mentum realizations, the conclusions listed above hold 
true. 
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